On the 1/c Expansion of f(R) Gravity 



O 

(N 

O 
O 



Joachim NaQ and Philippe Jetzer 
Institut fur Theoretische Physik, Universitat Zurich, 
Winterthurerstrasse 190, CH-8057 Zurich, Switzerland. 
(Dated: April 9, 2010) 

We derive for applications to isolated systems - on the scale of the Solar System - the first 
relativistic terms in the 1/c expansion of the space time metric <? M „ for metric f(R) gravity theories, 
where / is assumed to be analytic at R — 0. For our purpose it suffices to take into account up to 
quadratic terms in the expansion of f(R), thus we can approximate f(R) = R + aR 2 with a positive 
dimensional parameter a. In the non-relativistic limit, we get an additional Yukawa correction with 
coupling strength G/3 and Compton wave length y/6a to the Newtonian potential, which is a known 
result in the literature. As an application, we derive to the same order the correction to the geodetic 
precession of a gyroscope in a gravitational field and the precession of binary pulsars. The result 
of the Gravity Probe B experiment yields the limit a < 5 x 10 11 m 2 , whereas for the pulsar B in 
■ the PSR J0737-3039 system we get a bound which is about 10 4 times larger. On the other hand 

| the Eot-Wash experiment provides the best laboratory bound a < 10 _10 m 2 . Although the former 

bounds from geodesic precession are much larger than the laboratory ones, they are still meaningful 
in the case some type of chameleon effect is present and thus the effective values could be different 
' at different length scales. 

^ , PACS numbers: 04.25.-g; 04.25.Nx; 04.50. Kd 
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^ 1 I. INTRODUCTION 

Since the emergence of the concepts of dark matter (DM) and dark energy (DE), they still lack in a concrete and 
t-H ■ satisfying physical model. This open question motivated the development of new gravity theories. Most of them 
are direct modifications of general relativity (GR), which is still the simplest relativistic gravity theory fitting very 
accurately many precision measurements in astrophysics, such as Mercury perihelion shift or mass diagrams of double 
pulsars. Among such modified theories a lot of attention has been devoted to the so-called metric f(R) theories with 
an action 



-gd 4 x + S M , (1) 



where in contrast to GR the Einstein- Hilbert Lagrangian density is replaced by a nonlinear function f(R). Sm is the 
standard matter action. For an overview one may consult e. g. [l| and references therein. 

In the literature there are several approaches which address the question of the non-relativistic limit as well as 
relativistic approximations of metric /(i?) theories. A discussion of the first relativistic corrections after a transfor- 
mation to the Einstein frame is given in 0, Q • The non-relativistic limit in the Jordan frame is investigated in (343 , 
whereas a calculation in the Palatini formalism is given for example in [7J. In the present paper we work strictly in 
the Jordan frame. Our work is mainly motivated by the fact that the parametricpost-Newtonian (PPN) formalism is 
not adapted to cover the 1/c expansion of f(R) gravity @. As pointed out in @, 0|, the corresponding non-relativistic 
limit indeed is not Newtonian, but contains a Yukawa type correction, too. We therefore derive the lowest order 
relativistic terms of the 1 / c-expansion of the space time metric governed by the Euler Lagrange equations of ([TJ . We 
thus achieve a "post- Yukawa" approximation of f(R) gravity. This approximation is analogous to the complete first 
post-Newtonian approximation of GR, cf. for example [9l-[r2|. 

In section II we present the field equations of the model. Section III is devoted to the calculation of the expansion 
coefficients, and in section IV we make some remarks on the non-relativistic limit as well as the GR limit of the model. 
In section V we derive the equations of motion for a test particle and determine the underlying potentials for a set 
of freely falling particles. In section VI we apply our results to the precession of orbiting gyroscopes, and by using 
the measurements of Gravity Probe B and of the pulsar B in the PSR J0737-3039 system, we get upper limits for the 
value of a. 
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As far as notation is concerned: Greek letters denote space time indices and range from to 3, whereas Latin letters 
denote space indices and range from 1 to 3. We take the sum over repeated indices within a term. By an index ", /x" 
we denote the partial differentiation w.r.t to x 11 , except for /i = 0, where it denotes the differentiation w.r.t the time 
coordinate t rather than the coordinate a; = ct. 



II. THE FIELD EQUATIONS 

Consider a 4-dimensional Lorentz manifold with metric g^ v of signature (—,+,+,+). We write g — detg^ and 
denote the Ricci tensor of g M „ by The variation of the action ([lj w. r. t. the metric yields the Euler-Lagrange 

equations 

f(R)R»v - \}{R)9^ - V^ v f{R) + g^Dgf'iR) = ^-T^, (2) 

where R — g^R^v, — (—2c/ \/—g)(SSM / '5g^ v ) is the energy momentum tensor, c the vacuum speed of light, G 
Newton's constant, V M the covariant derivative for g^ and D g = V^V^. Taking the trace of ([2]) we obtain 

3D g f\R) + f'(R)R-2f(R) = ^-T, (3) 

where T is the trace of T^. Motivated by the post- Newtonian approximation of GR we calculate the coefficients of 
the expansion of g^ in powers of c _1 : 

ffoo = -l + (2) / lO o + (4) froo + 0(c- 6 ), (4) 

got = {3) h 0l + O(c- 5 ), 

9ij = <% + (2) ^-+0(c- 4 ), 

where denotes a quantity of order 0(c~ n ). The Ricci scalar is at least of order 0(c -2 ). Thus, if we assume 

the function / to be analytic at R = with /'(0) — 1, it suffices to consider the expansion 

f(R) = -2A + R + aR 2 , a ^ 0, (5) 

in order to calculate the coefficients of <? M „ up to the orders indicated in (j4]), since higher powers of R would only 
contribute to higher orders in the equations of the perturbation expansion. Moreover, since we adopt an expansion 
about a flat background space time in Q, we ignore a possible cosmological constant A in what follows. The influence 
of a nonvanishing A on the applications in section VI is discussed in [13]. As we will see later, the parameter a has 
to be positive for many reasons. 

It is convenient to introduce the scalar field <j) := f'(R)- Since f"(R) ^ holds for our choice of f{R) in (JSJ, we 
can invert f'{R) in order to show that f(R) gravity is equivalent to the Brans-Dicke theory with a non vanishing 
potential term and Brans-Dicke parameter lubd = 0. Wc define the scalar field tp by <p = 1 + 2atp, where we have 
chosen the asymptotic value such that a renormalization of the Newton's constant will be redundant (cf. the end of 
section IVA). Then the equations ^ and Q are equivalent to 

1 /8ttG / 1 \ 1 A 2 \\ 

Rfiv = 1 + 2aip I ~ ) + q9^^ + a I ^g^V + 2V M V M (^ I I (6) 

4ttG 1 

The field tp thus has the effective mass H/(cy/6a). From ([5]) we infer that the dimensionless quantity aR should be 
small compared to 1. This fact reflects the concept of the chameleon effect [T3|, which states the possibility that 
the Compton wave length A = \/6a of the field <p is smaller or larger in regions with higher or lower matter density, 
respectively. We understand our 1/c expansion to be valid in a local region which has an approximately constant 
mean matter density, in the sense that we assume the parameter a to be constant on the length scale characteristic 
for later applications, in particular the geodetic precession. On the other hand, a may vary for applications which 
have different length scales. 
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III. THE EXPANSION COEFFICIENTS 



We introduce space time coordinates (x, t), where bold face letters denote three dimensional vectors. The expansion 
coefficients are functions of these coordinates. Denote the 3-dimensional Nabla operator by V = (81,82, 83). We find 
the following expansions for and tp, 

Roo = ~ ( - V 2 Whoa - V 2 (4 ^oo + - {3) h 0l M - 4 (2) ^,oo (8) 



^ - V 2 ^hoo - V 2 Wftoo + - c i3) h 0h0l ~ ^ (2) h, 





+ (2) 


Roi — 


i( 


Rij = 


i( 


f = 


( 2 V 



U V 2 ^h tJ - Wh 00 ,ij ~ {2) h kkM + W hikM + <*) hkjtki + o ( c - 4 ), 



and for the energy momentum tensor, 

T 00 = (-2) T 00 + (0) T 00 + O ( c -2- )) (9) 

T oi = (-i) T oi + o( c - 1 ), 
T ij = (o) T ij + o(c~ 2 ). 

Equation ([7]) then yields in leading order the Yukawa equation 

v 2 ( 2 V- a 2 ( 2 V = -^^- 2 )T 00 . (10) 

c 4 

where we have defined a 2 := l/(6a) for a real a. Equation (TIU]) has the solution 



< a V(x,t) = 4nx,i) (11) 



-ex X— X 



2Ga 2 f (- 2 )T 00 (x',i)e^ ■ 



with the potential 

nx, t):=~-^-J - — | dV. (12) 

The 00-component of equation ([6]) at order C(c -2 ) is given by 

V M^ 00 = _^C- 2 ) T °° + :^V (13) 



and has the solution 



with the potentials 



< 2 >Mx, t) = i (2t/(x, i) - W(x, t)) (14) 



4G /■ (- 2 )T 00 (x',i) , 



U(*,t) :=^J |x _-^ V, (15) 
V(x',t) j3 , 



v ' ; 12tt 7 |x-x'| 

In contrast to iterated Coulomb integrals, the potential W is well defined because of the exponential decay of the 
Yukawa term in V . Before we determine the coefficients of the Oi- and ij-components, we impose the four gauge 
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conditions 

9ij,j ~ ~ (9jj ~ 9oo). t ~ = O (c~ 4 ) , (16) 

Using the condition (|16[) . we find for the ij-component of equation ([5]), up to order 0(c~ 2 ), 

V2(2) ^ = -(^? ( " 2)T °° + I (2 V)^ (18) 



with the solution 



< a >MM) = % (C/(x,t) + W(x,i)) . (19) 



Taking into account the gauge condition (JT7J), the Oi-component of equation ©, up to order 0(c 3 ), simplifies to 



Defining the potentials 



V 2 Whoi = ^^T« + ^ 2 )% , - ^ J /',,,,- (20) 



X (x,i) := ^ / l -- 1 r'" , (x'./)|x-x'|,/V. ,21) 



^(x,i) := -±- J y(x',t)|x-x'|dV, 



we can write equation (j2"U)) as 



V 2 ( (3) ^ + ^(x + ^ 0< )=^ ( - 2) T 0i . (22) 



The solution of equation (|2"2"j) is given by 

{3) h Ql ( x ,t) = i (y<(x,t) + Zi(x,t)) 

where 

4G ^(-DT 0i (x' J t) J 3 / 



(23) 



y<(x,t) := / — hi^dV, (24) 

c J |x-x'| 

Zi(x,t) := -(x(x,t)+V(x,t)) j0i . 

Thus, the /(-R) correction to the shift is only due to the gradient field Z. 

For the derivation of the component ^hoo, we first address ^tp. In view of the gauge (fT())) the 0(c -4 ) part of 
equation |J7J| is given by 

V 2 < 4 V - « 2 (4 V = (- (0)T ° + (25) 



+i ( Voo + (2 VV,« + ^ (a V 2 V, 



We use the identity 



U,iVi = i (w 2 {UV) - VV 2 U - C/V 2 U) (26) 
for two arbitrary potentials U and V as well as the equations (fT0|) . (jlll) and (fl9l) to rewrite equation (|25j) as 
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where 



A := 8nGa 2 ^^T 00 + ^T" - 1 (u + W - ^T 00 ^ (28) 



+V Q0 + a 2 UV + a 2 VW - —V 2 . 

From equation (|27[) we obtain 

1 

c 

with the potential 



<"V(x,t) = ^B(x,i), (29) 



^--s/ ^jgiP- ^+igyfr'). (30) 

The c _4 -component of equation ([6]) can be simplified with equation (jlTp : 

V 2 ^ o = <W« + (a) ftoo,i ( {2) h ijtj - \ 2 h jjti - ^h Q0 , ) (31) 



/ 2 (0) T 00 + (0) r » + ( (2) hm _ J_(2 )(fi ] (-2) T 00 



3c 4 V V 3a 

Using equations (dUJ), (HU), O, (HU, CGI]), (0, ^ and (gU) we write equation (JHU) in the form 

V (<^ + 1 (^ + J^V - > - £vw)) = 1 (32, 

with 



c = WnG f 2(0)r00 + (0)Tii + J_i 15U _ *y \ (_ 2 ) T oo j (33) 



3 V 2c2 V 3 " 2 

12 6 9cr 3 



Hence, 



where 



Wfcoo(x,t) = ll>(x,t), 



+7^(x,t)W(x,t) + ^ 3 V(x,i)W r (x,t). 



(34) 



t) (35) 



4 v ' ' v ' ' 6a 2 

The metric field, up to the hrst relativistic approximation, is thus determined by the fields U, W, Y, Z, D. 

IV. THE GR LIMIT AND THE NON-RELATIVISTIC LIMIT 
A. The GR Limit 

By taking the limit a — > rcsp. a — > oo, the theory converges to GR. Notice that with the assumption a > we 
have the following representation of the Dirac delta function: 

hm f I £(x') e , = / e(x')<5(x - x')dV = £(x) (36) 

q->oo Air x — x' / 
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for an arbitrary test function £ (x) . Hence 

Inn W(x, t) = g / ^.^^W = itf(x, t) (37) 

a->oo OC 2 J |X — X'| 2 

and thus 

lim (2) /i 00 (x,i) = lim (2) /in(x,t) (38) 

a— >oo a— »oo 

3 2G f (- 2 )T°y,t) , 3 , 2 

- ^(«.*) = — y [x _ x1 d x - ? M*,t), 

as expected. E/jv is exactly the Newtonian potential. There is no need for a rescaling of Newton's constant G, as it is 
sometimes the case for other modified gravity theories. 

B. The Non-Relativistic Limit 

From the equations (fT2")l . (fT4")l and (fT5|) we gather that for the f(R) model given by the non-relativistic limit is 
not Newtonian, since the component ^ /loo contains a Yukawa type term. 

More explicitly, if for instance we consider a perfect, non viscous fluid with mass density p, pressure p and velocity 
field v = (v 1 ^ 2 ^ 3 ), we have 

(-2)2*0 = ^2^ (3Q) 

("Dt 01 = pcv\ 

In the non-relativistic limit, the energy- momentum conservation 

T MV . V = (40) 

then yields the equation of continuity, 

d t p + d i {pv i ) = 0, (41) 

and the analogous of the Euler equation, 



p {d t v l + vtdjv*) = -d tP + pd, (u - ^W^j . 



(42) 



Aside from a Newtonian term, the potential W contains also a Yukawa type term. The gravitational force, represented 
by the second term on the right hand side, therefore contains also the gradient of a Yukawa potential. Together with 
the Poisson equations 

V*U = -^p, (43) 



(|4"Tj) and (|4"2"|) are the basic equations of the non-relativistic limit of hydrodynamics for metric f(R) theory. The 
parameter a = l/(6a 2 ) can be constrained by experiments which test a Yukawa type correction to the Newtonian 
potential. Experimental data and overviews can be found for example in [l5l - [20l |. Constraints on f(R) theories are 
given for instance in [21]]. From equation (|42l (see also equation (|49|) ) we find for our specific model the Yukawa field 
strength G/3. The Eot-Wash experiment [l5j thus yields the limit a < 10 _10 m 2 . 

V. PARTICLE DYNAMICS 

We can derive the equations of motion for a freely falling test particle in a field (U, W,Y, Z, D) by evaluating the 
geodesic equation, 



d 2 x^ dx v dx 



A 



Kx~j--j- = ^ (44) 



dr 2 ' ~ vX dr dr 
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where the denote the connection coefficients of the metric g^, and x^(t) — (ci(r),x(r)) is the position of the 
test particle at proper time r. Defining v := dx/dt, the equations of motion in vector notation read 

= v(u-lw)+^(U + W)v(u-\w)-VD (45) 

+<9 t (Y + Z) + v A (V A Y) - vd t Uu + ^1 

-3v(v ■ V)t7 - v 2 V(£7 + W)^j + O (c" 4 ). 
The energy momentum tensor of a set of point particles reads flp| 

"(x,i) = ^£m n ^^(^) V(x-x„(t)), (46) 



T^ 1 , 

dt dt \ dt 

n 

where m„ is the mass of the n-th particle, r„ is its proper time and x„(t) = (ir* (t), x 2 (i), x^(t)) is its position at time 
t. This leads to the expansion 



r 



on 



$^m„ (c 2 + ±(U~5W + v 2 n U8 3 ( X - Xn ) + 0(c- 2 ), (47) 

n ^ ' 

r« = £ ™««^ 3 (x - x„) + o (c- 2 ) , 

n 

where v n = cbc n /dt. Inserting (|47|) into the defining equations of the potentials U, W, x, V 7 an d Y{, we get 

X(x,t) = — ^m n |x-x„(t)| , 

n 

*(*,*) = |E^(|x-x4t)l + Q2|x _ 2 Xrt(t J l-e-^) 
„ |x-x„(*)| 

Using (j47]) and (J48]), one is able to calculate (at least formally) the potentials Z and D, which we however do not 
need for the applications we consider in this paper. 

We investigate the non-relativistic limit by taking the 0(1) part of equation (|451) . Considering each particle as a 
test particle in the field of the other ones, we replace the potentials U and W by their self-energy free parts. The 
equations of motion for the test particle then read 

This is the analogue of the Newtonian equations of motion for a purely gravitating set of point particles. 



VI. PRECESSION OF ORBITING GYROSCOPES 



The following derivation of the precession and its applications is done in complete analogy with the corresponding 
computations in GR fToL [TTj . A recent and detailed review of gravitomagnetism in physics and astrophysics is provided 
in 0. 
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Consider a gyroscope with spin S = (Si, 5*2, 5*3). We define its 3-velocity v to be the velocity of the non-relativistic 
center of mass of the gyroscope, the trajectory of which is assumed to be the one of a purely gravitating point particle. 
Notice that this last assumption includes the fact that the gyroscope moves along a geodesic and thus the Thomas 
precession due to a external force vanishes. The spin 4- vector := (Sq,S) precesses according to the equation of 
parallel transport, 



and satisfies the orthogonality condition 



= 0. (51) 



Up to the lowest order, equation (|50|) reads 

dS,, 



^ = (c (3) 4- (2) r^ + (2) r^)^. (52) 

Similarly as in GR we define the intrinsic spin vector by 

5:=(l-^(C/ + W))s-^v(vS). (53) 
Then S 2 is an integral of (1521 up to the required order. Equations (f5"2"j) and (|4"5j) then yield 

|^A5, (54) 



and the precession angular velocity is given by 



° = ~i (VAY) + ^ VAV ( C/+ 4 W '^ ^ 



Compared to GR, the Lense-Thirring precession represented by the first term remains unchanged. This was expected 
since, as mentioned in section III, a finite parameter a affects the shift only by a gradient field. On the other hand, 
the geodetic precession given by the second term is modified. 



A. Gyroscope Orbiting Around the Earth 



We now analyse the correction to the geodetic precession, since the Lense-Thirring precession is not modified. We 
model the Earth as a sphere with mass M which is at rest centred at the origin of our coordinate system. Consider 
the gyroscope to be in a circular orbit x(i) with radius |x(t)| = r and unit normal n, such that (x, v, n) is a positively 
oriented dreibein. Then (|48|) gives 

"+>=^fl-*H. (56) 



4 2r V 9 

and by equating the gravitational and centrifugal force on the gyroscope, we find for the velocity 

v = (t ( 1 + J( 1+ar ) e " r )) ; nAx - ( 57 ) 

Hence the geodetic precession angular velocity is 



^geodesic := jvAV ( U + -W) (58) 



-i-v A V (u + 7I 

3(GM) 3 / 2 / 1 

2c 2 r 5/2 \ + 3 [ 



— ar 



1/2 / 1 



1 - - (1 + ar)e- ar ) n. 
9 
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Obviously, 12 gco d C sic converges to its GR value for a — > oo. This result can be compared with the measurements of the 
Gravity Probe B experiment [23| . The measured value lies within a minimal residue of 30 mas/yr from the predicted 
GR value 6606 mas/yr. This allows to constrain the relative deviation from the GR value in (|58|) by approximately 
0.45%. Since this deviation decays faster than e~ ar with growing a, while r«7x 10 6 m, we expect a much larger 
bound for a as the one given by the Eot-Wash experiment. For the given accuracy of measurement, equation (|58[) 
yields a < 5 x 10 11 m 2 . If we estimate an upper limit for the scalar curvature using the mean Earth mass density, we 
are left with R < 10 -22 m~ 2 . Even in this very rough approximation we have aR <C 1 for our constraint on a. 

We remark that for a nonvanishing cosmological constant the residue in the Gravity Probe B measurements lead to 
the bound A < 3 x 10~ 27 m -2 [l3T |. When considering both A ^ and the Yukawa term within perturbation theory, 
to leading order the corrections due to a and to A will add linearly. 



B. Precession of Binary Pulsars 



Consider a binary system with center of mass at the origin. We index the mass m n and the position x„ by 1 for 
the pulsar and by 2 for the companion. Equation (|48|) gives for the fields of the companion 



^0 + >(M)- ^^(l-i^Ml), (59, 

Define now x := xi — X2, r := |x|, the reduced mass fj, := m\mij [m\ + 1712) and the angular momentum L := 
fix A (dx/dt). Evaluating (|55|) at xi then leads to 

n -S( 2+ £H (i+Qr H> (61 » 

In order to approximate the average of 12 over a period, (12) = T _1 J Q Sl(t)dt, we need an expression of the trajectory 
in the non-relativistic limit. Therefore we deduce the correction to the Kepler ellipse 

ro(0) := ,J fl , (62) 
1 + e cos 8 

due to a small perturbation SU of the Newtonian potential f/jy. In (|6"2"j) we have introduced polar coordinates (r, 9) 
and the parameters p :— L 2 /(Gmim 2 /i) and e := (1 + 2Ep j (Gm\ra-i)) x ^ > where £7 is the total energy of the two 
body system. We expand the equation of motion with respect to SU, 

^ ^^+^™+^V^V /2 (63) 



dff ILI 





-1 


)- 


r 2 y 


L 2 


v 1/2 




) 



/ Gm 1 m 2 \ l?\ x/ * / / Gm l m 2 \ L 2 \ 1/2 p t 
+0(SU 2 ). 

Integration along the unperturbed trajectory ro(0) for SU — Gmim2e~ ar / (3r) yields r(0) ps 7*o(0) + <5r((9) with 
We approximate the average of 12 as 

Due to the exponential decay the correction terms for r and (12) are very sensitive to variations of the parameter a. 
To give an idea of the orders of magnitude, we analyse (|6"5")) for the PSR 1913+16 data given in [24j. The predicted 
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GR value given by the first term on the r.h.s of (|65|) evaluates to 1.21°/yr. The correction term reaches approximately 
1% of the GR value for a sw 2.6 x 10 15 m 2 . Simultaneously, the correction 6r already after one period reaches 1% of 
the semi major axis, whose uncertainty can be measured with a much better accuracy of 4 x 10~ 5 %. From this last 
value we incidentally find the rough limit a < 1.7 x 10 14 m 2 , if we cumulate Sr for one year. 

The discovery of the double-pulsar binary PSR J0737-3039 paved the way to significantly improve the accu racy 
of binary pulsar measurements. An overview of the observed and derived parameters can be found in pM [26]. 
With this data we evaluate the precession rate for the pulsar B predicted by GR to 5.07°/yr. The measured value 
£Ib ~ 4.77^0 65 °/y r [13 then allows the correction to lie within a minimal residue of 7% from 5.07°/yr. This roughly 
yields the constraint a < 2.3 x 10 15 m 2 . 

For other binary pulsars, even the required accuracy to precisely test GR is not yet reached by the corresponding 
experimental research, see e. g. [28[ for the case of PSR J1141-6545. We finally conclude that a huge improvement of 
the accuracy of measurement would be necessary to put useful limits on a by the precession of binary pulsars. 



VII. CONCLUSIONS 



We gave the general formula for the lowest relativistic order coefficients of the 1 J c expansion for the metric of 
f(R) gravity, where we considered functions of the form f(R) = R + aR 2 . Furthermore, we investigated the GR and 
non-relativistic limits. The latter results in the Newtonian potential plus a Yukawa type correction with strength G/3 
and Compton wave length y/6a. As an application, we derived the /(-R)-corrections to geodetic precession of orbiting 
gyroscopes. The Lense-Thirring precession is not affected. 

While the laboratory bound from the Edt-Wash experiment provides the small bound a < 10~ 10 m 2 , the results 
from Gravity Probe B imply the much larger limit a < 5 x 10 11 m 2 . The measurements of the precession of the pulsar 
B in the PSR J0737-3039 system provide instead the limit a < 2.3 x 10 15 m 2 . Even for these large values of a the 
quadratic term in ([5]) still induces a small correction of GR. 

In principle, the coefficients for g^ v can be used for the same applications as the PPN coefficients of metric gravity 
theories which have a Newtonian non-relativistic limit. However, the computation of the applications which require 
the fourth order coefficient ^hoo are challenging, because its formula is quite involved and contains up to threefold 
iterated integrals. Therefore, a numerical analysis would be necessary for generic applications. 

It would be interesting to take into account more general functions /. For instance one could extend the choice of 
/ to functions which are not necessarily analytic at R = 0, but at a nonvanishing point R — Rq- Formally, this would 
imply to replace / given in ([5]) by the more general function 

f{R) = -2A + ai R + a 2 R 2 , A,a u a 2 ^0. (66) 

While the possibility of a\ ^= 1 would not cause much trouble in the derivation of the 1/c expansion, the nonvanishing 
cosmological constant requires an expansion about a de Sitter or anti-dc Sitter background, thus leading to more 
complicated partial differential equations for the potentials. Nevertheless, the choice of / as in would be needed 
to study many f(R) models which are proposed in the literature. 
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